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The growth of snow crystals is dependent on the temperature and saturation of the environment. 
In the case of dendrites, Reiter’s local two-dimensional model provides a realistic approach to the 
study of dendrite growth. In this paper we obtain a new geometric rule that incorporates interface 
control, a basic mechanism of crystallization that is not taken into account in the original Reiter’s 
model. By defining two new variables, growth latency and growth direction, our improved model 
gives a realistic model not only for dendrite but also for plate forms. 
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I. INTRODUCTION 

Snowflake growth is a specific example of crystalliza¬ 
tion - how crystals grow and create complex structures. 
Because crystallization corresponds to a basic phase tran¬ 
sition in physics, and crystals make up the foundation 
of several major industries, studying snowflake growth 
helps gaining understanding of how molecules condense 
to form crystals. This fundamental knowledge may help 
fabricate novel types of crystalline materials [7]. 

Whilst current computer modelling methods can gen¬ 
erate snowflake images that successfully capture some ba¬ 
sic features of real snowflakes, certain fundamental fea¬ 
tures of snowflakes growth are not well understood. One 
of the key challenges has been that the snowflake growth 
models consist of a large set of PDEs, and as in many 
chaos theoretic problems, rigorous study is difficult. 

Snowflakes exhibit a rich combination of characteristic 
symmetry and complexity. The six fold symmetry is a 
result of the hexagonal structure of the ice crystal lattice, 
and the complexity comes from the random motion of 
individual snow crystals falling through the atmosphere: 


FIG. 1: Examples of real plate and dendrite snowflakes [9]. 
(a) Stellar dendrite (b) Stellar plate (c) Sectored plate. 

Scientific studies of snowflakes can be categorized into 
two main types. The first approach takes a macroscopic 
view by observing natural snowflakes in a variety of mor¬ 
phological environments characterized by temperature, 
pressure and vapour density (see [12-14] and references 
therein). The second type takes a microscopic view and 
investigates the basic physical mechanisms governing the 
growth of snowflakes (e.g. see [7]). While some aspects of 
snowflake growth (e.g., the crystal structure of ice) are 
well understood, many other aspects such as diffusion 
limited growth are at best understood at a qualitative 
level [7]. 


Another approach in which snowflake growth is numer¬ 
ically simulated to produce images with mathematical 
models derived from physical principles is through com¬ 
puter modelling (e.g. see [3, 4, 15-17, 19]). By comparing 
computer generated images with actual snowflakes, one 
can correlate the mathematical models and their parame¬ 
ters with physical conditions. While computer modelling 
can generate snowflake images that successfully capture 
some basic features of actual snowflakes, so far there has 
been only limited analysis of these computer models in 
the literature. 

In this paper we analyze snowflake growth simulated 
by the computer models so as to connect the microscopic 
and macroscopic views and to further our understanding 
of snowflake physics. The models that have been con¬ 
sidered in the past are in essence chaos theoretic models, 
which is why they successfully capture the real world phe¬ 
nomena, but prove to be notoriously difficult to analyze 
rigorously. In this paper we study Reiter’s popular model 
[17] using a combined approach of mathematical analysis 
and numerical simulation. 

After reviewing Reiter’s model in Section II, in Section 
III we divide a snowflake image into main branches and 
side branches and define two new variables (growth la¬ 
tency and growth direction) to characterize the growth 
patterns. In Section IV we derive a closed form solution 
of the main branch growth latency using a one dimen¬ 
sional linear model, and compare it with the simulation 
results using the hexagonal automata. Then, in Section 
V we discover a few interesting patterns of the growth 
latency and direction of side branches. On the basis of 
the analysis and the principle of surface free energy min¬ 
imization, in Section VI we enhance Reiter’s model and 
thus obtain realistic results both for dendrites and plate 
forms. We summarize our contributions and present a 
few future work directions in Section VII. 


II. AN OVERVIEW OF REITER’S MODEL 

Reiter’s model is a hexagonal automata which can be 
described as follows. Given a tessellation of the plane into 
hexagonal cells, each cell z has six nearest neighbours. 
We shall denote by s t (z) G M>o the state variable of cell 
z at time t which gives the amount of water stored in z. 
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Then, cells are divided into three types: 

Definition 1 A cell z is frozen if s t (z) > 1 (an F-cell) . 
If a cell is not frozen itself but at least one of the nearest 
neighbours is frozen, the cell is a boundary cell (a B-cell). 
A cell that is neither frozen nor boundary is called nonre- 
ceptive (an NR-cell). The union of frozen and boundary 
cells are called receptive cells (R-cells). 


Frozen cell (F-cell) 
Boundary sell (B-cell) 
Nonreceptive cell (NR-cell) 
Edge cell (E-cell) 


FIG. 2: Classification of cells. 

The initial condition in Reiter’s model is 

1 if zmO 
(5 if z^O 

where O is the origin cell, and /3 represents a fixed con¬ 
stant background vapour level. 

Definition 2 Define the following functions on a cell z: 
the amount of water that participates in diffusion u t (z); 
and the amount of water that does not participate v t (z). 
Hence, 




s t (z) = u t (z) + v t (z) ( 1 ) 

and we let v t (z) := s t (z ) if z is receptive, and v t (z) := 0 
if z is non-receptive. 

For 7 , a two fixed constants representing vapour addition 
and diffusion coefficients respectively, in Reiter’s model 
the state of a cell evolves as a function of the states of its 
nearest neighbours according to two local update rules 
that reflect the underlying mathematical models: 

• Constant addition. For any receptive cell 2 , 

v t(z) :=Vt(z)+ 7 (2) 

• Diffusion. For any cell 2 , 

u t (z) ■= u i(z) + 2 (W (z) - u t (*'))> (3) 

where we have used upper indices ± to denote new func¬ 
tions giving the state variable of a cell before and after a 
step is completed, and written uf(z) for the average of 
uf for the six nearest neighbours of cell z. 

The underlying physical principle of Eq. (3) is the dif¬ 
fusion equation 


where a is a constant. Indeed, Eq. (3) is the discrete ver¬ 
sion of Eq. (4) on the hexagonal lattice, and it states that 
a cell z retains (1 — a/ 2 ) fraction of uf(z ), uniformly dis¬ 
tributes the remaining to its six neighbours, and receives 
a/12 fraction from each neighbour. The total amount 
of uffz) would be conserved within the entire system, 
except that a real world simulation consists of a finite 
number of contiguous cells. The cells at the edge of the 
simulation setup are referred to as edge cells , in which 
one sets u(~(z) := /3. Thus, water is added to the system 
via the edge cells in the diffusion process. Combining the 
two intermediate variables, one obtains 

s t+1 (z) :=u+(z)+v+(z). (5) 

By varying the parameters a, ffi 7 in Reiter’s model one 
can generate certain geometric forms of snowflakes ob¬ 
served in nature: 



(a) p = 0.3, 7 = 0.0001 (b) p = 0.35, 7 = 0.001 



(c) p = 0.6, 7 = 0.01 (d) p = 0.9, 7 = 0.05 


FIG. 3: Images generated by Reiter’s model [17] for a = 1. 


III. GENERAL GEOMETRIC PROPERTIES 

In what follows, we give new descriptions of snowflake 
growth and analyze them with a combined approach of 
mathematical analysis and numerical simulation by con¬ 
sidering a coordinate system of cells as in Figure 4(a) 
below. A cell z is represented by its coordinate for 

i,j £ Z, with the origin O = ( 0 , 0 ). Since there is a six 
fold symmetry, we only focus on one twelfth of the cells, 
marked as dark dots, for which j > i > 0 : 

The images in Figure 3 show that a crystal consists 
of six main branches that grow along the lattice axes, 
and numerous side branches that grow from the main 
branches in a seemingly random manner. The main and 
side branches exhibit a rich combination of characteristic 
symmetry and complexity which we shells rudy through 
the rate of water accumulation of a cell z, defined by 


( 4 ) 


A s t (z) = s t+ i(z) - s t {z). 


du/dt = a\7 2 u, 
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+ 90 ° 



- 90 ° 


+ 30 ° 

- 30 ° 


• z d 


• Z s 


(b) 


FIG. 4: (a) Coordinate system of hexagonal cells. 

(b) Definition of growth directions in the coordinate system 


Definition 4 Denote by Zd a destination cell, and by z s 
a source cell. Then, the growth direction of cell Zd is de¬ 
noted by g(zd) and defined as the orientation of z s with 
respect to Zd, where the angle is with respect to the hori¬ 
zontal axis. The source-destination cell relationship shall 
be denoted by S(z d ) := z s . 

As shown in Figure 4(b), the angle is given relative to 
the horizontal direction in the coordinate system, and 
satisfies 


At any time the set of all the R-cells is connected. 
Moreover since a frozen cell is surrounded by receptive 
cells and does not accumulate water via diffusion, and 
since water flows from nonreceptive cells to boundary 
cells, one has that the rate of water accumulation A s t (z) 
of a cell satisfies the following general geometric proper¬ 
ties: 

(A) For an NR-cell z, one has 0 < s t (z) < /3 and 
Ast(z) < 0. Suppose that an NR-cell 2 is sur¬ 
rounded by R-cells and disconnected from the E- 
cells. If 7 > 0, there exists to > 0 such that 
s t0 (z) > 1 (i.e., a time to in which the cell becomes 
frozen); otherwise, the NR-cell will permanently re¬ 
main non-receptive and never become frozen. 

(B) For a B-cell z the quantity A s t (z) is the sum of 
7 and diffusion. If 7 > 0, there exists to > 0 such 
that s to (z) > 1 ; otherwise, if cell z is surrounded by 
a set of F-cells and disconnected from the E-cells, 
then as in (A), the cell will never become frozen: 

lim s t (z) < 1 . 

t—yoo 

(C) For an F-cell z, one has A St(z) = 7. 

(D) The state variable of a cell is s t (z) = (3 only for z 
an E-cell. 

Thus, for an NR-cell to become frozen, the cell goes 
through two stages of growth. First, the NR-cell loses 
vapour to other cells due to diffusion, as in (A). Subse¬ 
quently, it becomes a B-cell and accumulates water via 
diffusion and addition, as in (B), until it becomes frozen 
and sees no benefit of diffusion, as in (C). Becoming a 
B-cell is a critical event between the two stages. 

We focus on the second stage and define two new vari¬ 
ables to characterize growth patterns. 

Definition 3 The time to be frozen of a cell z is denoted 
by T(z) and defined by the condition St( z ){ z ) > 1, and 
St(z ) < 1 for t < T{z). Similarly, we define B{z) as 
the first time to be boundary. Finally, growth latency is 
denoted by L(z) and defined by L(z) := T{z) — B(z). 

A cell becomes a B-cell as one of its neighbouring cells 
has just become an F-cell, and thus it is useful to make 
the following definition in terms of redistribution of wa¬ 
ter: 


g(z d ) G {+30°, -30°, +90°, -90°, +150°, -150°}. 

Note that while the growth of Zd is traced back to a 
unique z s , a source cell may correspond to multiple des¬ 
tination cells. 


IV. GROWTH OF MAIN BRANCHES 

Consider cells (i,j) where i + j = K for a fixed K. 
These cells are all K sites away from the origin (0, 0) on 
the grid. The main branch growth pattern is such that 
T(0, K) < T(i,j) and 

T (f’f) - T (^j) for even K, 

T ( K f 1 ’ K f k ) > T ihi) for odd K. 

Along the lattice j-axis, one has g(0, j) = —90° for all 
j. Hence, the snowflake growth is fastest along a lattice 
axis, which represents a main branch, and is the slowest 
along the 30°-offset lattice axis. 

We next develop a model to calculate the growth la¬ 
tency L( 0 , j). As cell ( 0 , j) becomes frozen, cell ( 0 , j + 1 ) 
becomes a boundary cell. Hence the first time to be 
boundary H(0, j + 1) = T(0, j), and thus one can calcu¬ 
late T(0, j) as 
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T(0,j)=T(0,0) + ]TL(0,fc). ( 6 ) 

k =1 

In order to gain analytical understanding, we first 
study a one dimensional model. Consider a line of con¬ 
secutive cells zo, zi, ..., zjv, where Zjsr is the edge cell. 
Initially cell O is frozen. We focus on the growth period 
[B(k),T(k)\ in which cells zo, zi, ..., Zk-i are frozen and 
cell k grows from boundary to frozen. Since Eq. (3) de¬ 
scribes the diffusion dynamics of vapour being transferred 
from the edge cell to cell Zfc, and cell z^ accumulates water 
via addition Eq. (2), to derive an analytical solution, we 
make the following assumption which we justify shortly. 

Assumption 1 For t G [B(zk),T(zk)\, assume that in 
Eq. (3) one has 

<{Zi) = Ut(Zi), 

for fc +1 < i < N, for N as above and B(zk),T(zk) as in 
Definition 3. Therefore, the vapour distribution reaches 
a steady state, denoted as /i(i\k). 
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From Assumption 1, we can ignore the notations of d=, 
and reduce Eq. (3) to the linear equation 

v(i\k)~^ ~ l\k) + + l\k)j . (7) 

Moreover, with the boundary conditions fi{k\k) = 0 and 
fi(N\k) = /3, the vapour distribution can be written in a 
closed form as follows: 

V(i\ k ) = f°r i - k, ..., N, (8) 

which graphically represents a line that connects the two 
boundary condition points. 

We shall now explain why Assumption 1 is well moti¬ 
vated. Suppose that the steady state distribution Eq. (8) 
is already reached at t = B(z k ), this is, 



FIG. 5: Comparison of vapour accumulation simulations for 
the parameters a = 1, ft = 0.4, 7 — 0.001. 


SB(z k )(zi) = n(i\k-l) = A—E /3 for i = k,...,N. 

Then, s t (zi) evolves in the interval of (B(z k ),T(zk)\ in 
the following manner. For N k, one has 

s B ( Zk ){zk) = »{k\k - 1 ) = + ~ 0 . 

Thus it is reasonable to assume 


L(z k ) = T(z k ) - B{z k ) » 1 


because cell z k will take several simulation steps to reach 
sr(z k )(zk) > 1- Moreover, 


\s B (z k )(zi) - n(i\k)\ 


i ~ (k ~ 1) g (i~k) 
N-(k-l) P ( N-k) 


< 1 


for A" >> k. Thus, in each simulation step for time 
t G (B(zk),T(zk)\, the function s t (zi) only varies slightly 
and can be considered approximately constant. Hence, 
ut(Zi) = Ut(Zi). 

From Eq. (3) and Eq. (7), we may estimate uf(zk) by 


ut(z k ) ■■= 


a 1 

4 N-k 


p. 


Moreover, since u t ( Zk ) =0 it follows that we can further 
estimate A s t (zk) and L(zk) by 


(z k ) 

L(z k ) 


f s B(z k )( z k ) _ 1 w-fe+l^ 

A§ t(z k ) fjvh^ + T'" 


( 9 ) 

( 10 ) 


One may also model L(zk) = T(zk) — B(zk) as a func¬ 
tion of cell index of the cells. In the one dimensional 
model with N = 50, Figure 6 compares L{k) determined 
by the simulation, and L{zk) predicted by Eq. (10) as 
the snowflake grows from the origin O to the edge cell. 
For any &, one observes that L(zk) < L(zk). This phe¬ 
nomenon is expected, since by solving the above PDE’s 
one has that there exists a > 0 such that at any time 
instance t G [B(zk),T(zk)\, for i = &,..., A, one has 

fi(i\k) < s t (zi) and A s t (z k ) < As t (z k ). 

As a result, L(z k ) > L(z k )- 



FIG. 6 : Comparison of growth latency simulations for the 
parameters a = 1, ft = 0.4, 7 — 0.001. 


In the one dimensional model with A = 50, we may 
compare the vapour accumulation in every simulation 
step, as the simulation proceeds from the time when cell 
k = 25 just becomes boundary to the time when it be¬ 
comes frozen. Figure 5 compares As t (z k ) at cell z k de¬ 
termined by the simulation, and A St(zk) predicted by 
Eq. (9) for time t G [B(zk),T(zk)\- Initially s t (k) 
/i(i|fc), and A s t (k) A s t (k). After about 5 simulation 

steps, A St(k) drops to a flat plateau, which is approxi¬ 
mately equal to sAk ). At any time £, one observes that 
A s t (k) < A s t (k). 


Equation Eq. (10) predicts that L(zk) drops monotoni- 
cally with k. In simulation, we observe that in the begin¬ 
ning the cells grow from boundary to frozen very quickly, 
well before the steady state is reached. As a result, the 
steady state Assumption 1 does not hold in that time 
period. Figure 6 shows that L{z k ) first increases, then 
drops, and eventually matches the prediction L(z k )- 
Finally, we return to the two dimensional hexagonal 
cellular case. With a similar steady state assumption, we 
can reduce the PDE to a set of linear equations similar to 
Eq. (7). However, the geometric structure is much more 
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complex than the one dimensional case. As a result, it 
is difficult to derive a closed form formula of the vapour 
distribution similar to Eq. ( 8 ). 

Figure 7 below plots L(0,j) along a main branch. 
Comparison with Figure 6 indicates a similarity between 
the one dimensional and two dimensional cases in that 
L(z) increases as the snowflake grows from the origin. 
However, in the two dimensional case, we observe from 
simulations that L(0,10) = L(0,11) = ... = L(0,195). 
When the snowflake grows close to the edge cell, it expe¬ 
riences some edge effect in the simulation where L drops 
drastically. This indicates that somewhat surprisingly 
Ast( 0 , j) remains almost constant as the snowflake grows 
along the main branch. 



FIG. 7: T(0, j) — R(0, j) of cells (0, j) along a main branch 

for j = 1,2,..., in the two dimensional scenario. Here cell 
(0, 200) is an edge cell, and a = 1, /3 = 0.4, 7 = 0.001. 


V. GROWTH OF SIDE BRANCHES 

While the main branches of snowflakes represent clean 
six fold symmetry, the side branches exhibit characteris¬ 
tic features of chaotic dynamics: complexity and unpre¬ 
dictability. Reiter’s model is completely deterministic 
with no noise or randomness involved, and yet the re¬ 
sultant snowflake images are sensitive to the parameters 
a, fi, and 7 in a chaotic manner. Chaos may appear to 
be the antithesis of symmetry and structure. Our goal 
in this section is to discover growth patterns that emerge 
from seemingly chaotic dynamics. 

Definition 5 Starting from a cell zq on the j-axis main 
branch, the set of consecutive frozen cells in the i-axis 
direction are referred to as side branch from cell zq. We 
shall denote by ze(zq) the outmost cell or tip, by E(zq) 
the length of the side branch, and the side branch itself 
by$(z 0 ) := {z 0 ,...,z E (z 0 )}. 

In what follows, we study the growth latency of side 
branches. Figure 8 below plots the tips of the side 
branches that grow from the j-axis main branch using 
the parameters of the four images in Figure 3. 



FIG. 8 : Plots of the tips (thin curves) and envelope curves 
(thick curves) of the side branches from the j-axis main 
branch using the parameters of the four example images in 
Figure 3. Due to symmetry, we focus on one set of side 
branches that grow from the right side of lattice j-axis. Here, 
the black curve represents Figure 3(a), blue for Figure 3(b), 
red for Figure 3(c), and magenta for Figure 3(d). The axes 
here are the horizontal and vertical axes of the coordinate 
system. 


Due to the chaotic dynamics, the lengths of the side 
branches vary drastically with zo in a seemingly random 
manner. For image (a), most of the side branches are 
short and only a small number stand out. The opposite 
holds for image (d). The scenarios are in between for 
images (b) and (c). The length of the side branches is 
indicative of the growth latency. The long side branches 
represent the ones that grow fastest. In Figure 8 we con¬ 
nect the tips of the long side branches to form an enve¬ 
lope curve that represents the frontier of the side branch 
growth. The most interesting observation is that the en¬ 
velope curve can be closely approximated by a straight 
line for the most part. Recall that the growth latency of 
the main branch is a constant. Thus we infer that the 
growth latency of the long side branches is also constant. 
Denoting by Lm and Ls the growth latencies of the main 
and long side branches respectively, one has that 

L m _ sin 27r/3 -0 

L s ~ sin 6 * 1 ’ 

where 0 is the angle between the envelope curve straight 
line and the j-axis. As a specific example, for the ma¬ 
genta curve, the envelope curve of the long side branches 
grows almost as fast as the main branch, such that 
6^ ~ 7 t/3 and the resultant image, appearing in Figure 
3(d), is roughly a hexagon. 

We shall consider next the growth directions of the cells 
on side branches. Figure 9 below plots the trace of the 
growth direction g(z) (see Definition 4) as a snowflake 
develops in the simulation. The corresponding snowflake 
image is shown in Figure 3(b). When a cell z becomes 
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boundary, we mark the cell to indicate g(z ) using the leg¬ 
end labeled in the figure. If a cell never becomes bound¬ 
ary, no mark is made. All side branches grow from the j- 
axis main branch, starting in the direction parallel to the 
i-axis. Subsequently, a side branch may split into mul¬ 
tiple directions. Indeed, all six orientations have been 
observed and the dynamics appear chaotic as g(z ) ap¬ 
pears unpredictable. However, we do find an interesting 
pattern described below. 


Compared with the straight paths, the deviating paths 
do not grow very far, because they compete with other 
straight or deviating paths for vapour accumulation in 
diffusion. On the other hand, the competition with the 
deviating paths slows down or may even block the growth 
of a straight path. When a straight path is blocked, the 
straight path is a strict subset of the corresponding side 
branch. This scenario is illustrated in Figure 10, where 
three side branches are shown. 



180 

160 

140 

120 

100 


Side branch cluseter 


i —Straight path 


j-axis 




Main branch 
Side branch 

Side branch cluster 

Destination cell 


Straight path 

Source cell • Destination cell 

Deviating path 


FIG. 10: An arrow linking two cells indicates the 

source/destination relationship. 


FIG. 9: Trace of relative orientations of source cells with re¬ 
spective destination cells. A destination cell becomes bound¬ 
ary because a source cell, which is one of the neighbours of 
the destination cell, becomes frozen. Legend is as follows: 
magenta ■ : +30°, black * : —30°, green + : +90°, 
blue o : —90°, red • : +150°, cyan X : —150°, where the 
parameters are a = 1, j3 = 0.35, 7 = 0.001. Note that not 
all straight paths are labeled. The axis are as in Figure 4(a). 


Definition 6 A straight path from a cell zq on the j- 
axis main branch is the set of consecutive frozen cells in 
the i-axis direction satisfying zi -1 = S(zi). The number 
of consecutive cells satisfying Zi -1 = S{zi ) is the length 
F{zq). We then denote a straight path from a cell zq by 

h^o) ; = {Z0,Z 1 ,Z 2 ,...,Z F(zo) }. 

Comparison between Definition 5 and Definition 6 
shows that the paths are nested, i.e. T(zo) C $(zo), and 
hence the lengths satisfy F(zq) < E(zq). When a cell Zi -1 
on the straight path becomes frozen, it triggers not only 
Zi in the i-axis direction but also other neighbours to be¬ 
come boundary, resulting in growth in other directions, 
which we call deviating paths. The straight and deviating 
paths collectively form a side branch cluster. 

Definition 7 The set of frozen cells that can be traced 
back to a cell on the straight path from cell zq on the j- 
axis main branch is referred to as a side branch cluster 
and denoted by ©(^o). 

A side branch cluster is a visual notion of a collec¬ 
tion of side branches that appear to grow together. Fig¬ 
ure 9 shows several side branch clusters and the cells 
on the corresponding straight path marked with cyan x. 


The straight path of the middle side branch is blocked 
by a deviating path of the lower side branch, which grows 
into a sizeable side branch cluster. Through the above 
definitions one has that if there exists a cell 2 such that 
z £ ©(zo) and z 0 3>(zo), then the paths are nested 
T(^q) C Moreover, the straight path determines 

the length of the side branch cluster: 

Definition 8 Denote by D(z,zq) the distance between 
zo,z £ ©( 20 ), defined as the smallest number of sites on 
the lattice between z and zq. The length of &(zq) is 

D(z 0 ) := max D(z,Zq). 
z e ©Go) 

Through the above definition one can show that there 
are K cells z% £ @(^o) such that the distances satisfy 
D(zi , zq) = D{zq) for i = 1,..., K with K > 2. Further¬ 
more, there exists zi £ T(zo) for 1 < i < K, and thus 
D(zq) = F(zq), the length of T(zo) as in Definition 6. 


VI. AN ENHANCED REITER’S MODEL 

Plates and dendrites are two basic types of regular, 
symmetrical snowflakes. We observe that while the den¬ 
drite images in Figure 3(a) (b) generated by Reiter’s 
model resemble quite accurately the real snowflake in 
Figure 1(a), as seen in Figure 3(c) (d) and Figure l(b)(c), 
the plate images differ significantly. The plate images in 
Figure 3(c) (d) is in effect generated as a very leafy den¬ 
drite. One of the reasons that Reiter’s model is unable 
to generate plate images realistically is that the model 
only includes diffusion, thus not taking into account the 
effect of local geometry. 
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As described in [1], two basic types of mechanisms con¬ 
tribute to the solidification process of snowflakes: dif¬ 
fusion control and interface control. Diffusion control 
is a non geometric growth model, where snowflake sur¬ 
faces are everywhere rough due to diffusion instability, 
a characteristic result of chaotic dynamics. For exam¬ 
ple, if a plane snowflake surface develops a small bump, 
it will have more exposure into the surrounding vapour 
and grow faster than its immediate neighbourhood due 
to diffusion. Interface control is a geometric mechanism 
where snowflake growth only depends on local geometry, 
i.e., curvature related forces. In the small bump example, 
the surface molecules on the bump with positive curva¬ 
ture have fewer nearest neighbours than do those on a 
plane surface and are thus more likely to be removed, 
making the bump move back to the plane. Interface con¬ 
trol makes snowflake surfaces smooth and stable, and it 
is illustrated in Figure 11 below. 


Destabilizing force 
A (diffusion control) 


Vapor region 



Snow crystal region 


FIG. 11: Two competing forces of diffusion control and in¬ 
terface control that determine snowflake growth. 

In summary, snowflake growth is determined by the 
competition of the destabilizing force (diffusion control) 
and stabilizing force (interface control). In the absence 
of interface control, Reiter’s model is unable to simulate 
certain features of snowflake growth. 

The interface between the snowflake and vapour re¬ 
gions has potential energy, called surface free energy, due 
to the unfilled electron orbitals of the surface molecules. 
The surface free energy 7 ( 71 ) as a function of direction n, 
is determined by the internal structure of the snowflake, 
and in the case of a lattice plane, is proportional to lat¬ 
tice spacing in a given direction. Figure 12 below plots 
the surface free energy y(n) of a snowflake as a function 
of the direction n. 



FIG. 12: Surface free energy of snowflake as a function of 
direction and equilibrium crystal shape of snowflake derived 
from surface free energy plot with Wulff construction [1]. 

The equilibrium shape of the interface is the one that 
minimizes the total surface free energy for a given en¬ 
closed volume. Wulff construction (see [1]) can be used 


to derive the equilibrium crystal shape W 7 from the sur¬ 
face free energy plot y(n): 

W 7 := {r | r • n < y(n), Vn}. (12) 

Wulff construction states that the distances of the equi¬ 
librium crystal shape from the origin are proportional to 
their surface free energies per unit area. Figure 12 plots 
the equilibrium crystal shape of snowflake. Moreover, it 
shows that due to interface control, snowflake growth is 
the slowest along the lattice axes, and the fastest along 
the 30°-offset lattice axes. 

This can be explained intuitively. Snowflake grows 
by adding layers of molecules to the existing surfaces. 
The larger the spaces between parallel lattice planes, the 
faster the growth is in that direction. This effect is com¬ 
pletely opposite to the diffusion control we have studied 
in Section 4, where snowflake grows fastest along the lat¬ 
tice axes. This is an example of competition between 
diffusion control and interface control. 

We next propose a new geometric rule to incorporate 
interface control in Reiter’s model. The idea is that 
the surface free energy minimization forces the lattice 
points on an equilibrium crystal shape to possess the 
same amount of vapour so that the surface tends to con¬ 
verge to the equilibrium crystal shape as the snowflake 
grows. 

From Figure 12, we learn that the equilibrium crystal 
shape is a hexagon except for six narrow regions along 
the 30°-offset lattice axes where the transition from one 
edge of the hexagon to another edge is smoothed. The 
equilibrium crystal shape used in the new geometric rule 
is shown in Figure 13. 



FIG. 13: Equilibrium crystal shape used in the new geomet¬ 
ric rule. 

Figure 13 shows the equilibrium crystal shape used in 
the new geometric rule and the interface control neigh¬ 
bours of the cells. As an example, cells A, B , C, E, E, F 
are on the same equilibrium crystal shape. Cells F,B are 
the interface control neighbours of A, cells A , C are the 
interface control neighbours of E, etc. 

The new geometric rule is applied after Eq. (5): a new 
variable S t (z) is defined to represent the amount of water 
to be redistributed for cell z at time £, with initial value 
St(z) = 0 for all z. 























Definition 9 For a given cell zo, define two interface 
control neighbours Zq,Zq, which are two neighbouring 
cells of zo on the same equilibrium crystal shape. 

Define 5 ( 2 : 0 ) as the average of the water amounts in cell 
zq and its two interface control neighbours Zq, Zq, this is: 



For every boundary 2 : 0 , if neither of Zq, Zq are frozen, 
then adjust St(zo) as follows 

$s(z 0 ) = S 8 (zo) +£( 5 ( 2 : 0 ) - sf +1 (z 0 )), (14) 

£ s 0o) = ^s(^o) “t” ^(^(^o) — s t+ 1 (^o))> ( 15 ) 

$s(4) = S s(4)+ e (H z o) ~ si +1 (z^)), ( 16 ) 

where 6 G M>o determines the amount of interface con¬ 

trol. After S s (z) has been adjusted for all z according to 
Eq. (14)-(16), for every cell z set 

s t+i( z ) := s t+i( z ) 4" ( z ). (17) 

Recall that in the original Reiter’s model, once water 
is accumulated in a boundary cell, water stays perma¬ 
nently in that cell. The new function Eq. (17) forces 
water redistribution particularly among boundary cells 
to smoothen the snow vapour interface. Figure 14 below 
shows two snowflake images generated by the enhanced 
Reiter’s model with the new geometric rule. 


(a) e = 0.1 (b) e = 0.01, a = 1, 0 = 0.4, 7 = 0.001 

FIG. 14: Snowflake images generated by the enhanced Re¬ 
iter’s model with the new geometric rule. 

At e = 0.1, the image above resembles a plate ob¬ 
served in nature much more closely than the ones in Fig¬ 
ure 3. By reducing interface control with 6 = 0.01, the 
snowflake starts as a plate and later becomes a dendrite 
as diffusion control dominates interface control. 


VII. CONCLUSIONS AND FUTURE WORK 

In this paper we have analyzed the growth of snowflake 
images generated by a computer simulation model (Re¬ 
iter’s model [17]), and have proposed ways to improve 
the model. We have derived an analytical solution of the 
main branch growth latency and made numerical com¬ 
parison with simulation results. Subsequently we ob¬ 
served interesting patterns of side branches in terms of 


growth latency and direction. Finally, to enhance the 
model, we have introduced a new geometric rule that in¬ 
corporates interface control, a basic mechanism of the 
solidification process, which is not present in the original 
Reiter’s model. 

The present work has shed light into some interest¬ 
ing patterns that lead to further questions about crys¬ 
tal growth. On the main branch growth, one may ask 
why the growth latency is almost constant (Figure 7) and 
whether this phenomenon is unique to the hexagonal cells 
or applicable to other two dimensional lattices. Concern¬ 
ing the side branch growth, it was noted that some side 
branches grow much faster than their neighbours, and 
that with slightly different diffusion parameters the side 
branch growth latency could change drastically at the 
same position while the main branch growth latency re¬ 
mains virtually the same. The study in Section VI shows 
that this great sensitivity is attributable to diffusion in¬ 
stability - when the growth of cells in some direction gain 
initial advantage over their neighbours, the advantage 
continues to expand such that the growth in that direc¬ 
tion becomes even faster. It was noted in Section VI that 
diffusion instability is caused by competition among cells 
in diffusion, and thus the average number of contributing 
neighbours is a good indicator to explain diffusion insta¬ 
bility. Finally, the enhanced model described in Section 
VI can be used to explore the interplay of diffusion and 
interface control. For example, one may simulate growth 
in an environment where the diffusion and interface con¬ 
trol parameters vary with time so as to generate images 
similar to Figure 1 (b)(c). 

Recently Reiter’s model was used in the study of snow¬ 
fall retrieval algorithms (e.g. see [6, 18] and references 
therein), and it was suggested that other mechanisms of 
snowflake formation from ice crystals besides aggregation 
must be considered in snowfall retrieval algorithms. It is 
thus natural to ask whether the enhanced Reiter’s model 
constructed here may provide novel insights in this direc¬ 
tion, as well as when considering crystal growth dynamics 
as in [8]. Moreover, since cellular automata models have 
been considered for numerical computations of pattern 
formation in snow crystal growth, it would be interesting 
analyze the outcome of the implementation of the model 
presented here to the analysis done in [2, 5]. 

Finally, the effects of lattice anisotropy coupled to a 
diffusion process have been studied in [10] to understand 
phase diagrams associated to crystal growth. Since this 
approach seemed recently useful from different perspec¬ 
tives (e.g. see [11] and references therein), it would be 
interesting to study the enhanced model constructed here 
from the perspective of [10, 11]. 
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